% idea is to put all thee time series together
% to see individual examples of substorm effects
load c:\manoj\projects\plasma\HUA_PIU1.mat;
load c:\manoj\projects\ace\kak_h;
load C:\Manoj\projects\ace\Julia_W_new.mat w_climate Julia_W fday w julia_model_drift ;
jul_fday = fday;
load c:\manoj\projects\ace\OMNI_ELEC.mat
% load c:\manoj\projects\ace\substorm;
load c:\manoj\projects\substorm\selected_substorm selected_substorm non_selected_substorm;
fday_substorm = non_selected_substorm;
%load c:\manoj\projects\ace\rand_substorm;
load  c:\manoj\projects\longp\final_coefficients a_n b_n ;
a = a_n;
b = b_n;
phase_delay = 17;%minutes
st = 18/24;
en = 20/24;
kak_mean = nanmean(KAK_H_DAILY);
Julia_fday = floor(jul_fday(:,1));
w=(w_climate)*24.366*1e-3; % climate removed w in mV/m

mean_hua = nanmean(HUAPIU);
flag1 = 0;
flag2 = 0;
flag3 = 0;
flag4 = 0;
flag5 = 0;
flag6 = 0;

radius_epoch =1/24;

L = fday_substorm >jul_fday(1,1);
k = find(L ==1);

for i = k(1):length(fday_substorm),
    
    
     if fday_substorm(i)-floor(fday_substorm(i)) >= st & ...
             fday_substorm(i)-floor(fday_substorm(i)) <= en,
%         

%        L = ace_all(:,1)  > fday_substorm(i)-  ( radius_epoch + 2/24)  & ...
%                ace_all(:,1) < fday_substorm(i) +  ( radius_epoch + 1/24);
%                  
%       ace_time = ace_all(L,1)+(phase_delay)/(60*24);%to delay ace data 
%       ace_data = ace_all(L,2);
%       
%       if all(~isnan(ace_data)),
%       
% 
%       pred_julia =filter(b,a,ace_data);
%       L = ace_time >= (fday_substorm(i) + 1e-5 ) - radius_epoch & ...
%               ace_time <= (fday_substorm(i) + 1e-5 ) + radius_epoch;
%           
%       SEG_PRED_JULIA = pred_julia(L); 
%       subplot(615);
%       plot(datenum(2000,1,1)+ace_time(L),SEG_PRED_JULIA);
%       ylabel('EEF mV/m');
%       hold on;
%       nk = axis;
%      line([datenum(2000,1,1)+fday_substorm(i) fday_substorm(i)+datenum(2000,1,1)],nk(3:4),'color','k');
%      hold off;
%      datetick('x',15)
%       flag1 = 1;
%       end;
         
    
       L = ace_all(:,1) > fday_substorm(i)-radius_epoch & ...
               ace_all(:,1) < fday_substorm(i)+radius_epoch;
                 
       if sum(L) > 20,
           DATA_STORM = ace_all(L,2);
           ace_storm_fday = floor(min(ace_all(L,1)));

           fprintf('FDAY = %d NUMBER OF DATA = %d\n',ace_storm_fday,sum(L));     
          
            if sum(L) > floor(radius_epoch*2*1440/5),
               kkk = find(L ==1);
               L(kkk(1)) = 0;
            end;
           
           
            
           SEG_ACE = ace_all(L,2);
           SEG_IMF = ace_gsm_bz(L);
           
           subplot(511);
           plot(datenum(2000,1,1)+ace_all(L,1),SEG_ACE);
           ylabel(' IEF Ey mV/km');
            hold on;
      nk = axis;
     line([datenum(2000,1,1)+fday_substorm(i) fday_substorm(i)+datenum(2000,1,1)],nk(3:4),'color','k');
     hold off;
     datetick('x',15)
           subplot(512);
           plot(datenum(2000,1,1)+ace_all(L,1),SEG_IMF);
           ylabel(' IMF Bz nT');
           flag2 = 1;
            hold on;
      nk = axis;
     line([datenum(2000,1,1)+fday_substorm(i) fday_substorm(i)+datenum(2000,1,1)],nk(3:4),'color','k');
     hold off;
     datetick('x',15)
       
           end;
           
          
  if fday_substorm(i)-floor(fday_substorm(i)) > st & ...
            fday_substorm(i)-floor(fday_substorm(i)) < en,
        
       [fday_int,ii] = intersect(fday_kak,floor(fday_substorm(i)));
       if ~isempty(fday_int),
          fday_frac = fday_kak(ii)+magnetic_frac;
           if fday_frac(1) < fday_substorm(i)-radius_epoch & ...
               fday_frac(end) > fday_substorm(i)+radius_epoch,

           
           L = fday_frac > fday_substorm(i)-radius_epoch & ...
               fday_frac < fday_substorm(i)+radius_epoch;
           
           if sum(L) > floor(radius_epoch*2*1440),
               kkk = find(L ==1);
               L(kkk(1)) = 0;
           end;
           
           dummy = KAK_H_DAILY(ii,:)- kak_mean;
           SEG_KAK = dummy(L);
           subplot(515);
           plot(datenum(2000,1,1)+fday_frac(L),dummy(L));
           ylabel('KAK H nT');
            flag3=1;
             hold on;
      nk = axis;
     line([datenum(2000,1,1)+fday_substorm(i) fday_substorm(i)+datenum(2000,1,1)],nk(3:4),'color','k');
     hold off;
     datetick('x',15)
           end;
           
           end;
       end;
       
       
          if fday_substorm(i)-floor(fday_substorm(i)) > st & ...
            fday_substorm(i)-floor(fday_substorm(i)) < en,
        
       [fday_int,ii] = intersect(HUA_PIU_fday,floor(fday_substorm(i)));
       
       if ~isempty(fday_int),
          fday_frac = HUA_PIU_fday(ii)+magnetic_frac;
           if fday_frac(1) < fday_substorm(i)-radius_epoch & ...
               fday_frac(end) > fday_substorm(i)+radius_epoch,

           
           L = fday_frac > fday_substorm(i)-radius_epoch & ...
               fday_frac < fday_substorm(i)+radius_epoch;
           
           if sum(L) > floor(radius_epoch*2*1440),
               kkk = find(L ==1);
               L(kkk(1)) = 0;
           end;
           HUAPIU(ii,:)= HUAPIU(ii,:)-mean_hua;
           SEG_HUAPIU = HUAPIU(ii,L);
           subplot(513);
           plot(datenum(2000,1,1)+fday_frac(L),SEG_HUAPIU);
           ylabel('\Delta H nT');
            hold on;
      nk = axis;
     line([datenum(2000,1,1)+fday_substorm(i) fday_substorm(i)+datenum(2000,1,1)],nk(3:4),'color','k');
     hold off;
     datetick('x',15)
            flag4=1;
           end;
           
           end;
       end;


   if fday_substorm(i)-floor(fday_substorm(i)) > st & ...
            fday_substorm(i)-floor(fday_substorm(i)) < en,
        
       [fday_int,ii] = intersect(Julia_fday,floor(fday_substorm(i)));
       if ~isempty(fday_int),
           
           if fday_substorm(i)-radius_epoch >= fday(ii,Julia_W(ii).k(1)) &...
               fday_substorm(i)+radius_epoch <= fday(ii,Julia_W(ii).k(end)),

           
           L = jul_fday(ii,:) >= (fday_substorm(i) + 1e-5 ) - radius_epoch & ...
               jul_fday(ii,:) <= (fday_substorm(i) + 1e-5 ) + radius_epoch;
           %1e-5 is to give a small offset to avoid fday_substorm at
           %exact time as a JULIA recording (hope so)
           
           %w(ii,:) = w(ii,:) - nanmean(w); % remove the daily forcing
           
           SEG_JULI = w(ii,L);
           subplot(514);
           plot(datenum(2000,1,1)+jul_fday(ii,L),SEG_JULI);
           ylabel('EEF mV/m');
            hold on;
      nk = axis;
     line([datenum(2000,1,1)+fday_substorm(i) fday_substorm(i)+datenum(2000,1,1)],nk(3:4),'color','k');
     hold off;
     datetick('x',15)
            flag5=1;
           end;
       end;
 
   end;
  
    
    %plot the time series _f
    
  
 if all([flag5]),
   % saveas(gcf,[sprintf('fday_%7.2f',fday_substorm(i)) '.fig'],'fig');
   title(datestr(fday_substorm(i)+datenum(2000,1,1)));
   A = datevec(fday_substorm(i)+datenum(2000,1,1));
   fprintf('%d\n',dayofyear(A(1),A(2),A(3)));
   pause;
end;
 
 flag1 = 0;
flag2 = 0;
flag3 = 0;
flag4 = 0;
flag5 = 0;
flag6 = 0;

 
end;
end;         
% plot((-radius_epoch*1440) + 2.5 : 5 :  (radius_epoch*1440) - 2.5 ,nanmean( SEG_ACE),'b');
% hold on;
% plot((-radius_epoch*1440) + 2.5 : 5 :  (radius_epoch*1440) - 2.5 ,nanmean( SEG_PRED_JULIA),'r');

